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ABSTRACT 



Shear and bulk viscosity and thermal conductivity for argon, krypton, xenon, 
and methane and the binary mixtures argon+krypton and argon+methane 
were determined by equilibrium molecular dynamics with the Green-Kubo 
method. The fluids were modeled by spherical Lennard- Jones pair-potentials 
with parameters adjusted to experimental vapor liquid-equilibria data alone. 
Good agreement between the predictions from simulation and experimental 
data is found for shear viscosity and thermal conductivity of the pure fluids and 
binary mixtures. The simulation results for the bulk viscosity show only poor 
agreement with experimental data for most fluids, despite good agreement 
with other simulations data from the literature. This indicates that presently 
available experimental data for the bulk viscosity, a property which is difficult 
to measure, are inaccurate. 

KEYWORDS: Viscosity; Thermal conductivity; Green-Kubo; Lennard- Jones; 
Molecular Dynamics; Molecular Simulation. 
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1. Introduction 



Transport properties play an important role in many technical and natural 
processes. With the rapid increase of available computing power, molecular 
simulation in combination with molecular modeling is becoming an interesting 
option for describing transport properties in regions where experimental data 
are not available or difficult to obtain. The calculation of transport coefficients 
by molecular simulation can be achieved by non-equilibrium molecular dynam- 
ics (NEMD) or equilibrium molecular dynamics (EMD). In NEMD, transport 
coefficients are calculated as the ratio of a flux to an appropriate driving force, 
extrapolating to the limit of zero driving force [T]. In EMD transport coeffi- 
cients are often calculated by the Green-Kubo formalism [2|3j. The choice 
between EMD and NEMD is largely a matter of taste and inclination, see e.g. 
[HfSfG] . There are numerous contributions in the literature in which both meth- 
ods were applied for the calculation of shear viscosity [I|7f51[Tmil|12j . bulk 
viscosity [Tn|14|15j . and thermal conductivity [7]I1U|11II16P!7] with comparable 
performance. To our knowledge the most comprehensive study on transport 
coefficients of the spherical Lennard- Jones fluid is reported in [TSfTO] . Despite 
of the large number of publications on simulation of transport properties with 
the spherical Lennard- Jones potential, not much effort seems to have been 
spent on a comparison of simulation results with experimental data of real 
fluids. Exceptions are the works of Michels et al. ^20j and Heyes et al. [TOfTT] . 
Michels et al. [20] compared the self-diffusion coefficient to the Chapman- 
Enskog theory and experimental data for krypton and methane at high den- 
sities, but thermodynamic properties were not considered. Heyes et al. PUIITT] 
simulated both transport and thermodynamic properties of argon+krypton, 
argon+methane, and methane+nitrogen, but properties of pure fluids were 
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not considered. In other cases pT][^ more complex molecules like ethylene, 
carbon dioxide, phenol, alkanes, or carbon tetrachloride were modeled by the 
spherical Lennard- Jones potential, which surely is an oversimplification. Simu- 
lation data on transport properties were compared there to experimental data, 
but not the thermodynamic properties. 

The interest here lies in the quantitative evaluation of the performance of 
Lennard- Jones models p3|, which have been optimized for the accurate pre- 
diction of thermodynamic properties, in the description of transport proper- 
ties. In this work EMD is used to carry out a comprehensive comparison, of 
shear and bulk viscosity and thermal conductivity of pure fiuids and binary 
mixtures of noble gases and methane. The molecular models are taken from 
[23] . These models were adjusted only to vapor-liquid equilibria and yield ac- 
curate descriptions of static thermodynamic properties over a wide range of 
temperatures and densities. Furthermore, they were recently applied for the 
prediction of diffusion coefficients of pure and binary mixtures of simple fluids 
over a wide range of temperatures and densities with good results [21]. 
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2. Theoretical background 



Transport coefficients are associated to irreversible processes, however, it is 
possible to describe irreversible processes in terms of reversible microscopic 
fluctuations, through the fluctuation dissipation theory [25j. In that theory, 
it is shown that transport coeflicients can be calculated as integrals of time- 
correlation functions of appropriate quantities There are different meth- 
ods to relate transport coefficients to time-correlation functions; a good review 
can be found in [26j. 

2.1. Shear and bulk viscosity 

The shear viscosity rjs, as defined in Newton's "law" of viscosity, describes the 
resistance of a fluid to shear forces. It refers to the resistance of an infinitesimal 
volume element to shear at constant volume [27]. The shear viscosity can also 
be related to momentum transport under the infiuence of velocity gradients. 
From a microscopic point of view, the shear viscosity can be calculated by 
integration of the time-autocorrelation function of the off diagonal elements 
of the stress tensor J^y fM2^ 

where V is the volume, ks is the Boltzmann constant, T the temperature, and 
< ... > denotes the ensemble average. The statistics of the ensemble average 
in Eq. ([T]) can be improved using all three independent off diagonal elements 
of the stress tensor J^y, J^^ and J^^. For a pure fluid, the component J^^ of 
the microscopic stress tensor Jp is given by 
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^ ^ ^ du(r -) 

jr = E-^<-f-EE-^-^- (2) 

i=l i=l j>i ^ ij 

Here i and j are the indices of the particles and the upper indices x and y 
denote the vector components of the particle velocities Vi. Eqs. ([T]) and ^ 
can be applied directly to mixtures. 

On the other hand, the bulk viscosity r^^ refers to the resistance to dilatation 
of an infinitesimal volume element at constant shape [2^. The bulk viscosity 
in polyatomic molecules, is related to a characteristic time required for the 
transfer of energy from the translational to the internal degrees of freedom 
[5U] . Moreover, the bulk viscosity plays an important role to describe ultrasonic 
wave absorption and dispersion [31]. From a microscopic point of view, the 
bulk viscosity can be calculated by integration of the time-autocorrelation 
function of the diagonal elements of the stress tensor and an additional term 
that involves the product of pressure p and volume V that does not occur 
in the shear viscosity, cf. Eq.([T]). In the NVE ensemble the bulk viscosity is 
given by [28)129] 

rib = f dt {{r/{t) - pVit)) ■ ( J-(0) - pVm). (3) 

The component J^^ of the microscopic stress tensor Jp is given by 

^ ^ ^ du(r- ) 

jr = E-^-^r-EE-s-^- (4) 

1=1 i=l j>i ^ tj 

The statistics of the ensemble average in Eq. (jl]) can be improved using all 
three independent diagonal elements of the stress tensor J^^, J™, J^^, and 
their permutations. In the case of mixtures Eqs. ([3]) and (HI) can be directly 
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applied. This equation is used in many simulation studies on the bulk viscosity 
e.g. [lOmiTMS] . 



2.3. Thermal conductivity 

The thermal conductivity A, as defined in Fourier's "law" of heat conduction, 
characterizes the capability of a substance for molecular transport of energy 
driven by temperature gradients. It can be calculated by integration of the 
time-autocorrelation function of the elements of the microscopic heat fiow , 
and is given by [28f29] 



VkeT'^ Jo 



oo 



dt (J^it) ■ J^(0)). (5) 



Here, the heat fiow 3q for a pure fiuid is given by 



1 ^ ^ ^ r du(ri ) 



(6) 



=1 j>i 



where Vj is the velocity vector of particle i and r^j is the distance vector 
between particles i and j. The term in squared parenthesis denotes the dif- 
ference between a dyadic product and the unitary tensor I multiplied by the 
intramolecular potential energy u{rij). This description of the heat fiow is not 
sufficient for binary and multi-component mixtures. In mixtures both diffusion 
and energy transport occur coupled [32], so that energy can be transported on 
a molecular level by diffusion or by heat transport. In a binary mixture with 
the components a and (3 the heat flow is given by p^|28] 



7 



J, = i E E (-f?v? - E E E E K' ^ ^ - 1 ■ • V? 

^fc=Q.i=l k=al=ai=lj>i ^ ij 

-E^'Ev^ (7) 

k=a i=l 

where h'^ denotes the partial molar enthalpy of component k. The computation 
of the heat flow in a binary mixture can, in principle, be accomplished in one 
simulation, however, here two separate simulations were preferred. One NpT 
simulation was performed for the computation of the partial molar enthalpies, 
corresponding to the enthalpic part of the energy flow, and another NVE 
simulation for the calculation of the autocorrelation function of the heat flow. 

2.4- Molecular models 

In this work, noble gases and methane are considered. These molecules ex- 
hibit rather simple intermolecular interaction so that the description of the 
molecular interactions by the Lennard- Jones 12-6 (LJ) potential is sufficient 
and physically meaningful for many technically relevant applications [.33j. The 
LJ potential u is defined by 




where a is the LJ size parameter, e the LJ energy parameter and the in- 
termolecular distance between particles i and j. The parameters a and e are 
taken from [23] and given in Table [H They were adjusted to experimental 
pure substance vapor-liquid equilibrium data alone. For modeling mixtures, 
parameters for the unlike interactions are needed. Following previous work of 
our group [2l|34|35j . they are given by a modified Lorentz-Berthelot combi- 
nation rule with 
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((Til + 0-22) 



(9) 



(712 = 



2 



and 



ei2 — ^ ■ veiled. 



(10) 



where ^ is an adjustable binary interaction parameter. This parameter allows 
an accurate description of the binary mixture data and was determined in 
previous work by an adjustment to one experimental bubble point The 
binary interaction parameters used in the present work are listed in Table [2l 

2.5. Simulation details 

The molecular simulations were performed in a cubic box of volume V con- 
taining A^=500 or A^=864 particles modeled by the LJ potential. The cut-off 
radius was set to = 5a, however, for very high densities where V is small, Vc 
was set to half of the box length, standard techniques for periodic boundary 
conditions and long-range corrections were used The simulations were 
started with the particles in a face-centered-cubic lattice with randomly as- 
signed velocities, the total momentum of the system was set to zero and New- 
ton's equations of motion were solved with a Gear predictor-corrector of fifth 



order [37]. The time step for this algorithm was set to At - yei/mi/cri = 0.001. 
All transport coefficients were calculated in the NVE ensemble, using equa- 
tions dH), ([3]), ([5]) and (171). The relative fluctuations in the total energy in 
the NVE ensemble were less than 10~^ for the longest run. The simulations 
were initiated in a NVT ensemble until equilibrium at the desired density and 
temperature was reached. Between 100 000 and 200 000 time steps were used 
for the equilibration depending on the state point. Once the equilibrium is 
reached, the thermostat was turned off and then the NVE ensemble invoked 
to calculate the transport coefficients by averaging the appropriate autocorre- 
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lation function. The length of the production period depended on density and 
temperature of the state point. At least 3 000 independent autocorrelation 
functions were used in the calculation of each coefficient of viscosity and 4 000 
in the calculation of each coefficient of thermal conductivity. In theory as Eqs. 
dl]), © and show, the value of the transport coefficient are determined by 
an infinite time integral. In fact, however, the integral is evaluated based on 
the length of the simulation. Therefore, the integration must be stopped at 
some finite time, ensuring that the contribution of the long-time tail [3H] is 
small. Figure 1 shows the behavior of the different autocorrelation functions 
and their integrals given by Eqs. ([T]), and for the most dense state 
points of argon for each transport property. As can be seen, all autocorrela- 
tion functions decay after 2 ps to less than 1 % of their normalized value. 
Later they oscillate around zero. To consider the effect of the long time tail, 
the calculation of the autocorrelations functions was extended to 5.4 ps for 
thermal conductivity and shear viscosity and to 6.5 ps for bulk viscosity. This 
was done because this autocorrelation function exhibits the largest fluctuation 
around zero attributable to long time correlation. The statistical uncertainty 
of the transport coefficients and thermodynamic properties were estimated 
using the Fincham's method [22] • For the calculation of the thermal conduc- 
tivity of mixtures it is necessary to include the partial molar enthalpies. For 
that purpose, Widom's test particle insertion [41] was taken using 2 000 test 
particles after each time step, 100 000 time steps for reaching equilibrium and 
300 000 for production. Our codes to calculate shear and bulk viscosity and 
thermal conductivity were successfully tested with the simulation results of 
Shoen et al. for viscosity ^12], Heyes [H] for bulk viscosity and Vogelsang et 
al. [T6] for thermal conductivity. 
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3. Results and Discussion 



In this section ttie prediction of shear and bulk viscosity and thermal conduc- 
tivity are compared pointwise with experimental data. For the shear viscosity 
the correlation of Rowley et al. [HIS], which is based on molecular simulation 
results, was also used. 

Figure 2 shows the results for the shear viscosity of argon, krypton, xenon 
and methane in comparison with experimental data. The data are reported at 
different temperatures and were taken from Vargaftik [12] for the noble gases 
and from Evers et al. [13] for methane. Overall, very good agreement between 
simulation and experimental data is found. The lowest relative deviations are 
found for argon and krypton with a few percent at lower densities. Also the 
results of shear viscosity of xenon and methane show very good agreement at 
low density, however, as the density increases, the deviations from the exper- 
imental data reach up to about 15% for xenon and about 20% for methane. 
It can be observed that the simulations for krypton, xenon and methane tend 
to underestimate the experimental viscosities as the density increases. This 
underestimation is larger in the results given by Rowley's correlation for the 
shear viscosity. 

Figure 3 shows the results for the shear viscosity of the binary mixtures ar- 
gon+krypton and argon-|-methane for two temperatures, the experimental 
data were taken from Mikhailenko et al. [i^flS] . Good agreement between sim- 
ulation and experimental data is found. For the mixture argon+krypton the 
typical deviations are about 10% and the highest deviations occur for krypton- 
rich state points. The predictions for the mixture argon+methane show a bet- 
ter agreement with the experimental data than those for argon+krypton. Typ- 
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ical deviations of the mixture argon+methane are about 10% at 100 K and 6% 
at 120 K, simulations performed at one intermediate temperature (T=110 K) 
confirm that better agreement between simulation and experiments is found 
as the temperature increases. The comparison of the present simulations with 
previous simulations of Heyes [TUfTT] . shows a comparable agreement for the 
mixture argon+krypton, however, in the mixture argon+methane the present 
simulations show better agreement than those of Heyes, specially in methane- 
rich state points, c.f. Fig. 3. 

Figure 4 shows the results for the bulk viscosity of argon, krypton, xenon 
and methane in comparison with experimental data. The experimental data 
are reported at different temperatures and were taken from Cowan et al, 
Malbrunotet et al, Cowan et al. and Singer [I6|47f48l|^ . respectively. The 
agreement is poor. Neither the density dependence nor the absolute value of 
rjiy predicted by molecular simulation agrees with the experimental data. The 
best results for rji, are achieved at the low temperatures and high densities for 
krypton. In this case the typical error is about 13%, even here the density 
dependence in not predicted correctly. For the other fluids, the predictions are 
lower than the experimental data by about 50%. Likewise, the experimental 
data of bulk viscosity show a stronger dependence on the density than the 
simulations. It must be pointed out here, however, that the method to measure 
the bulk viscosity by means of acoustic absorption of sound waves, involves 
considerable error |50J. Among the quoted experimental data, the krypton 
data are claimed to be the most accurate with an error band of about 25%, 
for the remainder error bands of up to 40% can be assumed. In the light of the 
fact that these simple molecular models describe both the thermal and caloric 
properties accurately [51] as well as the transport properties self-diffusion [21] , 
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shear viscosity and thermal conductivity (see below), it can be argued that 
the deviations founded for the bulk viscosity are due to experimental error. 

Figure 5 shows the results for bulk viscosity of the binary mixtures argons- 
krypton and argon+methane. The agreement is again poor. Neither the com- 
position dependence nor the absolute value of rjb predicted by molecular simu- 
lation agrees with the experimental data. The best results for rjb are achieved 
at the lowest temperature for the mixture argon-|-krypton. In this dis- 
crepancy of about 50% is found. In agreement with the results for pure fluids, 
it is found that the predictions are much too low in comparison to experimen- 
tal data. Previous work on these mixtures by Heynes et al. [TOlITT] confirms 
lower values from simulation, c.f. Fig. 5. 

Figure 6 shows the results for the thermal conductivity of argon, krypton, 
xenon and methane in comparison with experimental data. The data are re- 
ported at different temperatures along the bubble line and were taken from 
Vargaftik et al. [42J. Overall, very good agreement between simulation and ex- 
perimental data is present. The lowest relative deviations are found for argon 
and xenon, typical values are 4% for argon and 7% for xenon. The deviations 
for krypton and methane reach up to about 20%, no tendency is observed in 
these deviations. The good agreement for methane is especially remarkable, 
considering that the molecular model is very simplified and does not consider 
the contribution of rotation or internal degrees of freedom j52f53f54j . 

Figure 7 shows the results for the thermal conductivity of the binary mixtures 
argon-|-krypton and argon-(-methane for two temperatures, the experimental 
data were taken from Mikhailenko et al. [Hf^ . In general good agreement 
between simulation and experimental data is found. Due to the error intro- 
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duced by the partial molar enthalpy, the statistical uncertainty of the thermal 
conductivity was estimated as 5%. For the mixture argon+krypton the typical 
deviations are about 5% at 120 K, and about 7% at 140 K. For most of the 
simulated state points, these deviations lie within the uncertainty bars. For the 
mixture argon+methane a better agreement is found than for the mixture ar- 
gon+krypton. Over the whole composition range, simulation and experiment 
for argon+methane agree within the statistical uncertainties. The comparison 
of the present simulations with previous simulations of Heynes plOlITT] shows 
good agreement, c.f. Fig. 7. 
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4. Conclusion 



In the present work, the Green-Kubo formahsm was used to calculate transport 
properties for pure and binary mixtures of four noble gases and methane. The 
molecular interactions of the fluids were modeled by the spherical Lennard- 
Jones pair potential with parameters adjusted to vapor-liquid equilibrium only. 
A comprehensive comparison with available experimental data shows good 
agreement for pure fluids and binary mixtures for shear viscosity and thermal 
conductivity. On the other hand, for the bulk viscosity, with the exception 
of pure krypton, considerable systematic deviations between simulations and 
experiment occur. This disagreement hints towards highly inaccurate mea- 
surements. The present results support the finding that the spherical LJ 12-6 
potential is an adequate description for the regarded noble gases and also 
methane, in spite of the simplicity of the used model. Likewise the modified 
Lorentz-Berthelot combination rules with one binary interaction parameter 
are an adequate description of the molecular binary unlike interaction. It is 
worthwhile to extend the study to more complex fluids. 
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Greek Symbols 

a component 

(5 component 

At integration time step 

e Lennard- Jones energy parameter 

r]s shear viscosity 

r]y bulk viscosity 

A thermal conductivity 

^ adjustable binary interaction parameter 

a Lennard- Jones size parameter 

Vectorial and tensorial quantities 

I unitary matrix 

Jp stress tensor 

Jq heat flow vector 

Tij distance vector 

Vj velocity of particle i 
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Table 1 

Potential model parameters for the pure fluids used in this work [23] and molar 
mass [55J. 



Fluid 


a / k 


(e/A;B) / K 


M 1 g/mol 


neon 


2.8010 


33.921 


20.180 


argon 


3.3952 


116.79 


39.948 


krypton 


3.6274 


162.58 


83.8 


xenon 


3.9011 


227.55 


131.29 


methane 


3.7281 


148.55 


16.043 
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Table 2 

Binary interaction parameters taken from |34| . 



Mixture 




argon + krypton 


0.988 


argon + methane 


0.964 
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List of Figures 



1 Large plots: Autocorrelation functions. Small plots: Integrals 
following Eqs. ([T]), and All plots are shown for the 
most dense state points of argon for each transport property; 
top: thermal conductivity T=90 K and p=34433 mol ■ m~^, 
middle: shear viscosity T=150.7 K and p=35046 mol ■ m~^, 
bottom: bulk viscosity T=100 K and p=32843 mol ■ m~'^. 

2 Shear viscosity of argon, krypton, xenon and methane 
predicted by molecular simulation (full symbols) compared to 
experimental data (empty symbols) j^2l|l3] . argon T=300 K 
A; krypton T=230 K ■; xenon T=270 K methane 
r=100 K - 293.15 K correlation of Rowley et al. |8| - . 

3 Shear viscosity of the mixtures argon+krypton (top) and 
argon+metane (bottom) predicted by molecular simulation 
(full symbols) compared to experimental data (empty symbols) 
|4^E5] . The lines are a guide for the eye. 

4 Bulk viscosity of argon, krypton, xenon and methane 
predicted by molecular simulation (full symbols) compared to 
experimental data (empty symbols) [^47f48|im] . The lines 
are a guide for the eye. argon T=100 K - 145 K krypton 
T=116 K - 130 K ■; xenon T=165 K - 265 K methane 
T=100 K - 293.15 K ♦. 
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Bulk viscosity of the mixtures argon+krypton (top) and 
argon+metane (botton) predicted by molecular simulation 
(full symbols) compared to experimental data (empty symbols) 
(HE5] . The lines are a guide for the eye. 

Thermal conductivity of argon, krypton, xenon and methane 
predicted by molecular simulation (full symbols) compared to 
experimental data (empty symbols) [12]. The data correspond 
to bubble points reported at different temperatures, argon 
r=90 K - 140 K A; krypton T=140 K - 184 K ■; xenon 
T=170 K - 270 K methane T=100 K - 180 K ♦. 

Thermal conductivity of the mixtures argon+krypton (top) 
and argon+metane (botton) predicted by molecular simulation 
(full symbols) compared to experimental data (empty symbols) 
[HE5] . The lines are a guide for the eye. 
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Fig. 2. 
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Fig. 3. 




.0 0.2 0.4 0.6 0.8 1.0 



XAr / mol mol 



27 



Fig. 4. 
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Fig. 5. 
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Fig. 6. 
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Fig. 7. 
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